Expansion of Bose-Hubbard Mott insulators in optical lattices 
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We study the expansion of bosonic Mott insulators in the presence of an optical lattice after 
switching off a confining potential. We use the Gutzwiller mean-field approximation and consider 
two different setups. In the first one, the expansion is restricted to one direction. We show that 
this leads to the emergence of two condensates with well defined momenta, and argue that such a 
construct can be used to create atom lasers in optical lattices. In the second setup, we study Mott 
insulators that are allowed to expand in all directions in the lattice. In this case, a simple condensate 
is seen to develop within the mean-field approximation. However, its constituent bosons are found 
to populate many nonzero momentum modes. An analytic understanding of both phenomena in 
terms of the exact dispersion relation in the hard-core limit is presented. 
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I. INTRODUCTION 

In recent years, the study of strongly correlated ultra- 
cold gases in optical lattices has become an important 
focus of experimental and theoretical research Ini- 
tially, many of the efforts were devoted to studying equi- 
librium quantum phase transitions, for which the obser- 
vation of the superfluid to Mott-insulator phase transi- 
tions in three Q, two and one Q dimension was a 
major experimental accomplishment. More recently, the 
study of the noncquilibrium dynamics of these systems 
has started gaining attention [a-tlOj. The latter is possi- 
ble, and offers great insights into the coherent dynamics 
of quantum systems, thanks to the nearly ideal isola- 
tion of the gas from the environment, and to the unique 
control that has been achieved experimentally over the 
parameters that determine the dynamics. 

The noncquilibrium dynamics of quantum systems is 
a rich subject; see, e.g., Refs. [ITI - [l3| . In this paper, we 
will be mainly interested in the expansion of a strongly 
interacting bosonic gas in the presence of an optical lat- 
tice. Experimentally, expansions (time-of-flight measure- 
ments) have been used to learn about the properties of 
equilibrium gases 



141. The idea behind those measure- 



ments is that, if the gas is allowed to expand in the ab- 
sence of interactions, the initial momentum distribution 
function fully determines the density distribution after a 
long expansion time. Hence, by taking a picture of the 
latter, one can determine the former by means of a simple 
rescaling [l|, [l,[l4[- In such experiments, the particles are 
allowed to expand in absence of both the trapping poten- 
tial and the underlying optical lattice. This implies that 
the kinetic energy is significantly larger than the inter- 
action energy and, therefore, interactions can generally 
be neglected during the expansion. On the other hand, 
in the presence of the optical lattice where interactions 
cannot be neglected during the expansion, unusual phe- 
15| . e.g., the transformation of the 



nomena can occur 



momentum distribution of an expanding gas of impene- 
trable bosons into the momentum distribution function 



of a noninteracting Fermi gas in equilibrium, which has 
a Fermi edge [H], [13] • In addition, it has been recently 
proposed that the expansion of two-component fermions 
in the presence of strong interactions in an optical lattice 
can be used as a tool to cool the gas through a "quantum 
distillation" process fl8( | . 

There is another remarkable phenomenon that occurs 
during the expansion of a strongly interacting gas of 
bosons in an optical lattice, and that is the emergence of 
coherence during the expansion of Mott-insulating states. 
This phenomenon will be the main focus of the present 
paper. A Mott insulator is an interaction-generated in- 
sulator that exhibits a gap to one-particle excitations. 
The presence of this gap produces an exponential de- 
cay of the one-particle correlations. Within the Bose- 
Hubbard model description, the Mott insulator is the 
ground state of lattice bosons for commensurate fillings 
and sufficiently strong interactions (l9| . As mentioned 
before, Mott insulators have been created experimentally 
using ultracold gases in optical lattices. Due to the pres- 
ence of a trapping potential, the Mott insulating phases 
in experiments usually coexist with superfluid domains. 
Surprisingly, when Mott insulators are allowed to expand 
by turning off the trap in the presence of a lattice, co- 
herence emerges between initially uncorrelated particles 
[20M22I ] . a phenomenon that resembles a sort of dynami- 
cal phase transition. 

The study in Ref. [2(| demonstrated that the expan- 
sion of pure Mott-insulating states of hard-core bosons 
in one dimension (ID) leads to the development of quasi- 
long-range correlations and to the emergence of quasi- 
condensates at finite momenta k = ±7r/2a (where a is 
the lattice parameter). In a related setup, the onset of 
quasi-long-range correlations has been studied more re- 
cently (in the spin- 1/2 language) using various analytic 
and approximate approaches [231 ] . It is important to note 
that these results were obtained in the limit of infinite 
on-site repulsions between the bosons (hard-core limit), 
and for pure Mott-insulating domains. It is, therefore, 
natural to wonder whether the same results can be ob- 
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served experimentally with ultracold gases, for which the 
on-site interactions are not infinite and the insulating do- 
mains arc always surrounded by superfluid ones. Systems 
containing hard-core bosons, in which superfluid and in- 
sulating domains coexist, were studied in Ref. (2f[. It 
was found there that, if most of the initial system is in 
a Mott-insulating state, quasicondensatcs with momenta 
k ~ ±7r/2a are generated. The superfluid wings sur- 
rounding the insulating domain in the initial density pro- 
files appeared to have no negative effect on the emergence 
of quasicoherence during the expansion. Furthermore, 
the expansion of a Mott insulator in a system with fi- 
nite on-site interactions (soft-core bosons) was studied in 
Ref. [22|. Once again, quasi- long-range coherence devel- 
oped during the expansion and lead to quasicondensates 
with momenta |fc| < n/2a. The differences between the 
position of the peaks for soft-core bosons and hard-core 
bosons can be understood by considering their respective 
dispersion relations. Finally, a similar onset of power- 
law decaying correlations has been observed during the 
expansion of a Mott insulator within the fcrmionic Hub- 
bard model 0. 

A more challenging question is what happens in higher 
dimensions. This question poses greater conceptual and 
computational challenges. Since hard-core bosons in ID 
are described by an integrable model, one may expect 
their behavior to be unique. Hence it may not ex- 
tend to nonintcgrablc models far from any integrable 
limit, such as the Bosc-Hubbard model with finite on- 
site interaction in higher dimensions. In the particular 
case of one spatial dimension, since the boson's on-site 
interaction has to be strong enough for a Mott insu- 
lator to exist, the results obtained for the ID expan- 
sion [22| can be thought as still affected by the infinite- 
repulsion integrable limit, and thus may not generalize 
to higher dimensions. Furthermore, while the ID expan- 
sion can be studied by means of exact analytical and / or 
numerical approaches in the hard-core limit [2(| [23| 
and by means of an unbiased numerical technique for 
the soft-core case (the time-dependent-density-matrix- 
renormalization-group approach (25|), no such tool ex- 
ists to study the expansion in higher dimensions. For the 
latter, no efficient and unbiased, analytical, or computa- 
tional approach is known that can deal with the dynamics 
in the presence of strong interactions. Approximations 
are therefore required to study these systems. 

A first step toward understanding what happens in 
higher dimensions was taken in Ref. [26]. In that 
work, which was a study in the hard-core limit utiliz- 
ing the Gutzwillcr mean-field approximation, a three- 
dimensional (3D) Mott insulator was allowed to expand 
in a single direction in the lattice. The study showed 
that condensates emerge at finite momenta following the 
"melting" of a Mott insulator. In addition, the momenta 
of the resulting condensates were found to be fully de- 
termined by the ratio between the hopping amplitudes 
transverse to and along the expansion. Hence these 
systems can be used to create highly controllable atom 



lasers. Other noncquilibrium setups that allow one to 
obtain condensation at finite momentum [27], HH do not 
seem suitable for the latter purpose. 

In this paper, we extend the results reported in 
Ref. to the expansion of Mott insulators with fi- 
nite on-site interaction in two dimensions (2D). Our re- 
sults in 2D can be straightforwardly extended to three 
dimensions. We first consider a setup similar to that in 
Ref. [26[ , in which soft-core bosons are allowed to expand 
in a single direction. We show that similarly to the case 
in the hard-core limit, condensates emerge at finite mo- 
menta. We then proceed to study a system in which the 
Mott insulator is allowed to expand in all directions of the 
lattice. For this particular case, we show that, at least 
within the mean-field approximation, a simple conden- 
sate with many different momenta develops as the Mott 
insulator expands. We also discuss how these results can 
be understood in terms of the exact dispersion relation 
in the hard-core limit. 

The paper is organized as follows. In Sec. UH we intro- 
duce the model and the time-dependent mean-field ap- 
proach used to study the expansion. The explanation of 
how our initial states are prepared is presented in Scc. lIIII 
Section [IV] is devoted to the study of the expansion along 
one direction. The case in which the Mott insulator is 
allowed to expand in all directions is discussed in Sec. [V] 
Finally, the conclusions are presented in Sec. EH 

II. MODEL AND MEAN-FIELD APPROACH 

Our study is relevant to ultracold bosons trapped in 
deep optical lattices. They can be well described by the 
Bose-Hubbard Hamiltonian [29[ 

# = - £ J v ( 6 l a J + H - c -) + ^ - x ) 

+ Y^{V x x1 + V y yl)n i , (I) 

i 

where a\ (at) is the boson creation (annihilation) opera- 
tor at a given site i, and n, = a]a 4 is the particle num- 
ber operator. The first two terms of the expression are 
the usual kinetic and interaction terms that define the 
Hubbard model in the homogeneous case, where is 
the hopping amplitude between neighboring sites i and j 
and U is the on-site interaction strength. In ex- 
periments involving ultracold gases, a trap is required to 
confine the atoms. To a good approximation, those traps 
are harmonic, as considered in the last term in Eq. (TJJ. 
Here, V x and V y denote the strength of the trap in the 
x and y direction, while x% and yt are the coordinates 
of site i with respect to the center of the trap. In the 
remainder of the manuscript, positions will be given in 
units of the lattice spacing a, which is taken to be the 
same in all directions. 

The ground-state phase diagram of the homogeneous 
case (V x = V y = and Jy = J) has been studied exten- 
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sively in all dimensions (l9l. [30l - l35| . As mentioned earlier, 
it contains both superfluid and Mott-insulating phases. 
The former occurs for incommensurate fillings, and for 
commensurate fillings for U / J is below some critical value 
that depends on the filling and dimensionality of the sys- 
tem. The insulating phases are only present for commen- 
surate fillings and sufficiently strong interactions. This 
model is not exactly solvable in any dimension. Unbiased 
results for the phase diagram have been obtained using 

"Hill, density- 
, and strong cou- 
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quantum Monte Carlo simulations 
matrix rcnormalization group (ID) 

pling expansions (3lL |32| . Interestingly, a simple approx- 
imation such as the Gutzwiller mean-field theory, which 
is exact only in the limit of infinite dimensions, can still 
describe qualitatively the complete phase diagram of this 
model for finite dimensions [19j . 

In the presence of a trap, superfluid and Mott- 
insulating phases coexist and form a "wedding cake" 
structure [36l438j . The particular shape of the density 
distribution function can be understood within the local 
density approximation given the phase diagram in 
the homogeneous case (see Fig. [TJ and the accompanying 
explanation). Experimentally, such a structure has been 
observed by several groups [39|-|43| . As mentioned in the 
Introduction, the unique opportunities provided by the 
study of these systems come from the experimentalist's 
high degree of control over the depth and structure of 
the optical lattice and the confining potential. Essen- 
tially, all the microscopic parameters in Eq. ([1]) can be 
manipulated and made time dependent. For example, J, 
U, and V Xl V v can be modified by changing the inten- 
sity of the laser beams that produce the lattice; V x ,V y 
can be modified by introducing an additional trapping 
(or antitrapping) magnetic field; and U can be modified 
using Feshbach resonances. Together with the high de- 
gree of isolation from the environment, this level of con- 
trol makes these experiments ideal to study the effects of 
strong correlations and dimensionality in the nonequilib- 
rium dynamics of the expanding gases. 

In this paper, we are interested in studying the expan- 
sion of Mott-insulating states in the lattice. For this, we 
prepare the system in its ground state, for a particular 
number of particles, U, J, and V x ,V y . At time t = 0, 
we turn off the trap and allow the system to expand in 
the presence of the optical lattice. During the expansion, 
particles interact strongly and we compute observablcs 
such as the density, momentum distribution function, 
and condensate fraction, as a function of time. All our 
calculations are performed utilizing the time-dependent 
Gutzwiller mean-field theory. This is a good approxima- 
tion in high dimensions and, in Ref. |26j |. it was shown 
to qualitatively reproduce the behavior seen in an ex- 
act diagonalization study of the expansion in (small) 2D 
lattices. 

This mean-field theory is based on the restriction that 



the wave function is of the Gutzwillcr-typc product state, 



I* 



MF 



i=l n=0 



(2) 



where L is the number of lattice sites, n c is a cutoff taken 
to be large enough such that all our results are indepen- 
dent of n c (n c = 3 in our calculations), \n)i are the Fock 
states for lattice site i, and the complex coefficients ji n 
are variational parameters determined by minimizing the 
expectation value 



(^mf\H-^N\^ 



MF 



(3) 



In Eq. (|3"j). fx is the chemical potential, i.e., we work on 
the grand canonical ensemble, and N the total number 
of particle operator. 

The energy minimization leads to the set of equations 



-(n+l) + Vi-n 



+ n 



= 0, 



(4) 



where ^ ^ . denotes a summation over all j that are 
nearest neighbors of i, Vi = V x + V y yf, and 7i(-i) = 



7i(n c + l) 



0. The mean- field is defined as 



n=l 



<n lo{n _ l)l3n . 



(5) 



By solving Eq. (UJ, one obtains all the coefficients ji n 
in Eq. ([5]), and hence determines the initial state for the 
time evolution. 

An alternative approach, which is equivalent to that 
of solving Eq. (QJ, can be used. It consists in finding 
the ground state of the following mean-field decoupled 
Hamiltonian 1441: 



H 



MF 



^3 

(id) 



H.c. 



hi(fii - 1) + V" (Vi - /jl) hi. (6) 



This Hamiltonian is obtained by a decoupling of the hop- 
ping terms in Eq. ([!]) as 



a a,j 



(7) 



Equation (jTJ) can then be written as a sum over single-site 
Hamiltonians Hi, coupled only through constant neigh- 
boring terms $j. This implies that the ground state of 
the decoupled Hamiltonian can be written as the prod- 
uct state in Eq. ([2]). The optimal set of coefficients 7;„ 
is related to the ground-state eigenvector components of 
the Hamiltonian in Eq. ((SJ) subject to the condition in 
Eq. which has to be self-consistcntly satisfied. In 
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practice, this condition is reached as follows: an arbitrary 
set of nonzero $° is used to construct a set of decoupled 
Hamiltonians T-Li which arc diagonalized for each site in 
order to obtain a set of 7i„'s. Given j^ n , a set of $? ew is 
computed according to Eq. ([5]). If the elements on this 
new set are equal to those of the previous set within a 
certain desired precision, then self-consistency has been 
reached. If not, we set <£>° = $°™ and construct and a 
new set of decoupled Hamiltonians and the process of di- 
agonalizing and finding $^ ow is iteratively repeated until 
self-consistency is satisfied. We have verified that both 
approaches based on Eq. (j4]) and Eq. (j6]) provide quanti- 
tatively consistent results. However, at least within our 
implementations, the one based on the mean-field decou- 
pled Hamiltonian is more stable and considerably faster 
than the one based on the direct solution of Eq. (Ql . 

After turning off the trap, the time evolution of IV'mf) 
is computed using the time-dependent variational princi- 
ple [45l448| . which minimizes the expectation value, 

(*MF\id t - H + hN\*mf) , (8) 
and leads to the following set of differential equations: 
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FIG. 1. (Color online) (a) Mean- field phase diagram of the 
Bose-Hubbard model in 2D. The red arrow follows the dif- 
ferent phases occurring in the trapped system as one moves 
away from the center of the trap, (b) The density profile as 
obtained from the LDA with the largest Mott domain plotted 
together with the spatially varying chemical potential fXi/U. 



U 



(n-l) + Vi-n 



7m • (9) 



Here, we use the same notation as in Eq. (j4]). It is impor- 
tant to note that Eq. ([9]) determines the time evolution 
while preserving the normalization of the wave function 
and the total number of particles Nt in the system. We 
solve this set of L x (n c + 1) differential equations numer- 
ically, using a fourth-order Runge-Kutta method. Self- 
consistency is enforced by monitoring the total energy, 
particle number, and the normalization. 



the trap to be the largest chemical potential of the ho- 
mogeneous system that allows for a Mott insulator with 
density n = 1. That is just the value of the chemical po- 
tential fi + right at the transition from the Mott insulator 
with density re = 1 to the superfluid with n > 1. This 
is based on the observation that, as the distance from 
the center of the trap is increased, the effective chemical 
potential of the homogeneous system decreases. In 2D, 
the critical value ^ + can be determined within the mean- 
field theory by finding the largest solution for /i from the 
following relation that determines the boundaries of the 
Mott lobes [H, 



III. THE INITIAL STATE 

We are interested in the expansion of Mott-insulating 
states, though we know a priori that, in presence of a 
trapping potential, such Mott states are inevitably sur- 
rounded by superfluid regions. Nevertheless, one can aim 
at finding a setup for which those Mott domains are as 
large as possible for a given value of J/U. Our approach 
to achieving large Mott domains is based on the local- 
density approximation (LDA). The LDA assumes that 
the local properties of the confined system can be mapped 
to those properties of the homogeneous system with the 
same value of J/U and an effective chemical potential 
corresponding to the spatially varying /Xj = fx — Vi at site 
i in the trap. One can, therefore, get an approximate 
picture of what the trapped system looks like based on 
the phase diagram of the homogeneous system. Within 
the LDA, in order to get the largest Mott domain, it 
is enough to set the chemical potential at the center of 



2(J X + J V ) = jn-fi/U) (p/U-n+j) 
U 1 + fi/U 



where J x and J y are the hopping parameters along the x 
and y directions, respectively. These ideas are illustrated 
in Fig. [TJ where the mean-field phase diagram of the ho- 
mogeneous 2D system is depicted [Fig. (Ha)], as well as 
the density profile across the center of a trapped system 
that contains the largest Mott domain as obtained from 
the LDA [Fig. Hp))]. 

We should stress that, once fi + has been identified, we 
select a value of fj, close to /i + and solve Eq. Q and/or 
Eq. ^ to obtain the initial state. Hence our initial state 
is always the exact ground state within the mean-field ap- 
proximation, i.e., it is not the LDA ground state. This is 
the way the initial wave function in Eq. @ is constructed 
for all our expansions. 
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FIG. 2. (Color online) Comparison between the time evo- 
lution of the density profile across the trap of systems with 
U/J x = 25 (left) and U/J x = 60 (right), at t = 0, 10, 20, 30 
and 40 (from top to bottom) following the release from traps 
with Vx/Jx = 0.03 and V x /J x = 0.08, respectively. rj = 0.3 
and the time is reported in units of h/J x . 



IV. EXPANSION IN ONE DIRECTION 

Let us first consider the expansion of a Mott insula- 
tor along one direction in the lattice. We study the case 
in which the initial state is confined only along the x 
direction, i.e., V y = 0, while V x is kept finite, but the 
system has open boundary conditions in the y direction 
(a box trap). We allow for the hopping amplitudes to 
be anisotropic, namely, that the hopping along the x di- 
rection (J x ) differs from the one along the y direction 
(J y ). We then study the expansion for several values of 
U/J x , large enough so that the Mott-insulating phase ex- 
ists, as well as different values of rj = J y /J x . In all cases 
considered in this section, the expansion takes place in a 
lattice with 200 x 40 sites. As stated in the Introduction, 
this study extends the results presented in Ref. [26| on 
the expansion of hard-core bosons in three-dimensional 
lattices to the expansion of soft-core bosons confined in 
two-dimensional lattices. The initial state is selected such 
that the Mott-insulating domain is the largest possible 
for given values of U/J x , rj, and V x /J x , as explained in 
the previous section. 

In Fig. [2J we show the time evolution of the den- 
sity profiles of systems with U/J x = 25 and 60, for 
rj = 0.3. In the initial state (t = 0), both systems exhibit 
large Mott-insulating domains in the center of the lattice 
[x = 0), surrounded by small superfluid regions, the lat- 
ter being larger for the smallest interaction strength, i.e., 
U/J x = 25. For t > 0, the particles are allowed to ex- 
pand by switching off the confining potential along the x 
direction. As a result, the Mott-insulating domains melt. 
This process is seen to be qualitatively similar for the two 
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FIG. 3. (Color online) Comparison between the time evolu- 
tion of the momentum distribution function of systems with 
U/J x = 25 (left) and U/J x = 60 (right), at t = 0, 10, 20, 30 
and 40 (from top to bottom) following the release from traps 
with V x /J x — 0.03 and V x /J x = 0.08, respectively, rj — 0.3 
and the time is reported in units of h/J x . 

values of U / J x considered in Fig. [5] Note that during the 
time evolution, due to the symmetry of the initial state, 
the system expands only along the x direction, and the 
density profile maintains the reflection symmetry with 
respect to both x and y axes. 

The corresponding momentum distributions of the 
systems we have just described are presented in Fig. GO 
At t = 0, although barely observable in the plots (first 
row in Fig. GO), there is a sizable population of low momen- 
tum modes, because of the presence of the inevitable su- 
perfluid regions surrounding the initial Mott insulators. 
Those modes are embedded in a roughly homogeneous 
background of low-populated momenta. The latter are 
the result of the large Mott-insulating domains present 
in the initial states (top panels in Fig. [2]). Remarkably, 
shortly after the expansion has started, particles get re- 
distributed in momentum space, and they tend to fill 
up states surrounding two well-defined momentum values 
k = (±fc*, 0). This effect, which can be clearly observed 
in all panels with t > in Fig. GO is already present at 
very early times t < 10. As the expansion goes on, the 
values of the momentum of the highly populated modes 
become time independent. Their finite momenta reflects 
the motion of the two groups of particles moving away 
from the center of the systems in Fig. [2j We find that 
the structure around the momentum peaks, better seen 
in the left column in Fig. GO around k, disappears as the 
value of U/J x is increased, and only sharp peaks are left. 
This is apparent by comparing the momentum distribu- 
tion functions between t = 20 and t — 40 in Fig. GO There, 
the peaks for U/J x = 60 are much sharper than those for 
U/J x = 25. 
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To a very good approximation, on can predict the val- 
ues of k seen in Fig. [3J using energy and momentum 
conservation arguments. In order to do that, we first 
realize that, if U/J x is large enough, double or higher 
occupancies of the lattice sites become energetically too 
costly and can be, to a first approximation, ignored. The 
system can then be thought as composed of hard-core 
bosons. Following Ref. [26j, the kinetic energy can then 
be written as 



k 



£k"k, 



(11) 



where 7ik is the momentum distribution function at any 
given time, k = (k x , k y ), and £k is the dispersion relation 

(12) 



gk = — 2 J x cos k x 



*2tJy COS ky . 



In order to further simplify the calculations, one can 
assume that the system before the expansion is a pure 
Mott insulator, thus neglecting the small supcrfluid re- 
gions that are always present in the trap, and implying 
that the initial kinetic energy is -Ekin ~ 0. (Note that 
if double occupancies are neglected, the kinetic energy is 
essentially the total energy of the system after the trap 
has been turned off.) We further recall the observation 
that, in Fig. [3j many particles tend to populate modes 
with k = (±fc*,0), from which we understand that the 
y component of k is zero due to the symmetry of the 
expansion, which occurs only along the x direction. We 
can now ask, within our simplified picture of the systems 
presented in Figs. [2] and [31 which modes [with (fc*,0)] 
could be populated by all the particles after the expan- 
sion while still conserving the energy. (We are assuming 
that interactions during the expansion redistribute the 
energy in such a way that all particles condense into a 
single mode.) From Eq. (fTTj) . it follows that 



cos fc* 



-'/■ 



(13) 



If, in addition to energy conservation, we enforce con- 
servation of the total lattice momentum, or the symme- 
try of the expansion for that matter, we conclude that 
k = (±A;*,0), where k* = arccos(— rj). 

From Eq. ([T3"|). it follows that, in the limit of U/J x > 1, 
the location of the k modes that are populated during the 
expansion depends only on the ratio 77. Remarkably, we 
find that the peaks in Fig.[3Jagree very well with that pre- 
diction [Eq. ([T5]) ] . This occurs despite the fact that, in 
Fig. [3j the results were obtained for finite values of U / J x 
and in the presence of small supcrfluid domains around 
the Mott insulator. In Fig. Eta), we report a complete 
study of the dependence of the location of the peaks in 
the mean-field calculations as a function of 77, for three 
values of U/J x (points), together with the predictions of 
Eq. ([13)) (solid line). For the smallest values of U / J x that 
allow for the realization of a Mott insulator, one expects 
some dependence of the location of the peaks on U/J x , 
because the dispersion relation is expected to exhibit the 
largest differences from Eq. (|T2^) . The deviations of the 
mean-field results from Eq. (|13|) are, however, very small 
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FIG. 4. (Color online) (a) Location of the momentum peak 
k x after a long expansion time obtained from the mean- 
field calculation as a function of the hopping ratio 77, for 
U/J x = 25, 60, and 100 (points). The prediction of Eq. (TT3]) is 
depicted as a solid line, (b) The time evolution of the location 
of the momentum peak for U / J x = 60 and several values of 77. 
The horizontal dashed-dotted lines show the values predicted 
in Eq. ()13[) . The time is reported in units of h/J x . 



for small values 77. The largest deviations, for all values of 
U I J x , are observed for the largest values of rj. This occurs 
because as 77 is increased, the supcrfluid wings surround- 
ing the Mott insulator increase their size; as a result, the 
total kinetic energy is lower than that of the Mott insula- 
tor, and one of our assumptions (Ekin ~ 0) starts to fail. 
Nevertheless, as Fig. 0|a) shows, for the largest values 
of 77, the differences are still relatively small even for the 
smallest values of U/J x that support a Mott insulator. 
They also can be seen to decrease when increasing U / J x . 
For U/J x = 100, the results are barely distinguishable 
from the analytic results for the hard-core limit. Over- 
all, the fact that only for large values of U/J x one can 
realize a Mott insulator in the initial state allows for the 
hard-core description to be a good approximation in all 
cases analyzed here. 

We conclude this section by showing, in Fig. @|b), re- 
sults for the location of the momentum peak as function 
of time, for UJ J x = 60. That figure makes clear that, as 
time passes, the location of the peak rapidly approaches 
the values predicted by Eq. (|T3|) (dotted horizontal lines) . 
It also shows that it takes a shorter time to reach those 
values for smaller values of rj. In our calculations, the 
value of k* is determined as the local center of mass of 
the momentum distribution function around the position 
of the brightest peak. 



V. EXPANSION IN ALL DIRECTIONS 

We now examine the situation in which the initial Mott 
insulator is allowed to expand in all directions. For such 
a case, we consider V x — V y — V and J x = J y = J. We 
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FIG. 5. (Color online) Comparison between the time evo- 
lution of the density profile across the trap of systems with 
U/J = 25 (left) and U/J = 50 (right), at t = 0, 4, 8, and 
12 (from top to bottom) following the release from traps with 
V/J = 0.053 and V/J = 0.161, respectively. The time is 
reported in units of hj J. 
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FIG. 6. (Color online) Comparison between the time evolu- 
tion of the momentum distribution function of systems with 
U/J = 25 (left) and U/J = 50 (right), at t = 0, 4, 8, and 
12 (from top to bottom) following the release from traps with 
V/J = 0.053 and V/J = 0.161, respectively. The time is 
reported in units of hj J. 



turn off the trap at t = 0, and study the time evolution of 
the density and momentum distribution for systems with 
U/J = 25, 30, 40, 50, 60, 80, and 100. Note that the 
Mott insulator can only be achieved in systems with large 
enough values of U/J (U/J> 23.2 within the mean-field 
approach in 2D). 

Figure [5] depicts the time evolution of the density pro- 
file of the expanding bosons with U/J = 25 (left) and 
U / J = 50 (right). Interestingly, although the rightmost 
case has an interaction energy that is twice that of the 
leftmost one, the expansion is somehow similar in both 
cases. At the moment the trap is turned off (t = 0), the 
Mott-insulator domain is easily discernible from the su- 
pcrfluid one surrounding it. The lattice sites in the vicin- 
ity of the center of the system have unit occupancies and 
are surrounded by a thin ring with density smaller than 
one. Note that, as expected from the way we construct 
the initial state (see Sec. IIIip . the superfiuid ring in the 
system with U/J = 50 is smaller than the one in the 
system with U / J = 25. 

After turning off the trap, the bosonic cloud expands 
and the Mott insulator melts. In Fig. [51 it is remarkable 



that the expansion along the diagonals is always faster 
than that along the x and y axes. This leads to a square- 
like density profile for t > 0. Such a behavior is accen- 
tuated as U/J increases, as seen in the right panels in 
Fig. [5j where the expansion occurs almost entirely along 
the diagonals. 

The evolution of the momentum distribution function, 
corresponding to the density profiles depicted in Fig. [SJ 
is shown in Fig. [5] At t = 0, peaks can be seen in the 
center of the momentum distribution function. They cor- 
respond to the bosons that are in the superfiuid rings 
surrounding the Mott insulator. The difference between 
the size of those rings for U/J = 25 (Fig. [3J left) and 
U j J = 50 (Fig. [5j right) is apparent in the initial mo- 
mentum distribution, as the height of the peak for the 
U/J = 25 case is significantly larger than that in the 
U/J = 50 case. Surprisingly, for t > 0, these peaks 
rapidly evolve into circular-like structures, which in turn, 
become square-like as time passes (Fig. [5]). The initial 
ring structure is less evident in the U/J = 50 case at 
t = 4 because the formation of the square-like struc- 
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FIG. 7. (Color online) Top row: Comparison between the density profiles of systems with U/J = 25, 30, 50, 80, and fOO 
(from left to right), at t = 12 following the release from traps with V/J = 0.053, 0.0764, 0.161, 0.278 and 0.356, respectively. 
Bottom row: Corresponding momentum distribution functions. The time is reported in units of h/J. 



ture occurs faster as the value of U / J is increased. Also, 
the structure discernible in the momentum distribution 
function for small values of k x and k y (bottom panels 
for U/J — 25) disappears as the value of U/J is in- 
creased (see how much weaker it is in bottom panels for 
U/J = 50) and most of the bosons have momenta that 
lay within the edges of the square. The existence of this 
structure, in which bosons remain with small values of k x 
and k y , can be directly related to the size of the initial 
supcrfluid domain, which, as said before, decreases as the 
system approaches the hard-core limit. 

A comparison between the density and momentum 
distributions of systems with different values of U/J 
(U/J = 25, 30, 50, 80, and 100), at the time the bosons 
reach the boundaries of our 80 x 80 lattice (t = 12), 
is presented in Fig. [7] These plots make evident that 
the expansion mainly occurs along the diagonals as Uj J 
increases and, at the same time, the momentum distri- 
bution function acquires a clearer square shape. 

An understanding of the behavior exhibited by these 
systems, when allowed to expand in all directions, can 
be gained following a similar line of reasoning as done in 
the previous section. If U/J 1, a good description is 
provided by the hard-core boson picture. For hard-core 
bosons, after the trapping potential is turned off, the to- 
tal energy is the kinetic energy, just as in the expansion in 
one direction. The dispersion relation in the isotropic 
case reduces to 

£k = — 2 J(cos k x + cos k y ). (14) 

We now assume, once again to simplify the argument, 
that the system before the expansion is a perfect Mott in- 
sulator, i.e., all the particles are localized in their respec- 
tive sites, and occupy all possible momenta in k space. 
Therefore, the kinetic energy is again E^ ul = 0, and it 
is conserved during the expansion. Furthermore, let us 
assume that, after long time, the majority of bosons will 
occupy ( "condense" to) modes with particular k vectors 



(as seen in the mean- field calculations). Energy conser- 
vation is guaranteed by the assumption that the energy 
of each occupied mode vanishes. Under those assump- 
tions, the k vector of the modes populated by the bosons 
can be calculated from the dispersion relation as 

£k = — 2J(cosk x + cos k y ) = 0. (15) 

The solution to this equation is given by the square-like 
contour plot seen in Fig. [FJ for the case with w = 0. 
Analytically, such a contour plot is described by the ex- 
pressions 

k y = 7r ± k Xl ky = — 7r ± k x . (16) 

Remarkably, the momentum distribution of all systems 
in Fig. [7] has converged to a distinctive square-like struc- 
ture. It closely resembles the prediction of Eq. (fl"6|) . as 
represented in Fig. [8j The deviations from the exact re- 
sults in Eq. (fTfJ)) can be understood to be the result of (i) 
the finite values of U/J in all our calculations, and (ii) 
the fact that the Mott-insulating domains with n = 1 arc 
always surrounded by supcrfluid ones. Once again, the 
first modifies the dispersion relation in Eq. (fl"4")) , and the 
second one reduces the energy from the value assumed in 
Eq. (|T5|) . As clearly seen in Fig. when U/J increases 
and the supcrfluid domains in the initial state reduce 
their size, the momentum distribution becomes closer to 
the square-like contour defined in Eq. ([TBI . The robust- 
ness of the results, despite points (i) and (ii) above, is 
remarkable. Hence we expect that they should be re- 
producible in experiments with ultracold bosons even at 
finite, but low, temperatures. 

The analysis above also allows one to understand the 
expansion pattern seen in the density profiles in Figs. [5] 
and [71 The density profiles for t > arc determined 
by the velocities with which bosons expand across the 
lattice. In the hard-core limit, the group velocity is 

v 9 = Vek = 2 J(sin k x x + sin k y y), (17) 
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FIG. 8. (Color online) Depiction of the group velocity vector 
field, given by Eq. (|17[) in fc-space, along with a series of con- 
tour plots satisfying the equation — 2(cos k x + cos k v ) = w/J, 
where w/J = -2, -1.5, -1, -0.5, 0, 0.5, 1, and 2. 



which produces the vector field plot seen in Fig. [5] It is 
clear that the group velocities, which are perpendicular 
to the contour plots, increase in magnitude for those k 
modes along the diagonal. (The maximum group veloc- 
ity occurs when = along the diagonal.) Therefore, 
if such modes are highly populated, the net effect is that 
the system will expand preferably along the directions 
indicated by those group velocity vectors. Indeed, the 
high population of such modes is apparent in the mo- 
mentum distribution functions depicted in Fig. [7l The 
low population of the modes in the vertices of the square 
indicates the suppression of particles with momenta close 
to (±7r, 0) and (0, ±7r), which have vanishing group ve- 
locities. In Fig. [7J the accuracy of this description can 
be seen to improve with increasing U/J. It is interesting 
to note that the system expands at the greatest speed 
possible by populating the modes with the highest group 
velocity that are consistent with the conservation of the 
energy, as can be concluded from the momentum distri- 
butions and comparing them with the analytical results. 

We have already shown that, during the expansion 
from a mostly Mott-insulating state, bosons tend to pop- 
ulate momentum modes around well-defined values. Re- 
sults have been presented for mean-field calculations and 
supported by analytic arguments. One can now won- 
der whether all these highly populated modes are coher- 
ently coupled forming a simple Bose-Einstein condensate, 
or whether they are part of fragmented condensates, of 
quasicondensates, or their occupation is simply 0(N®) 
but their number increases with increasing system size. 
In strongly correlated systems, it is possible that modes 
with different momenta couple to each other and still 
form a simple Bose-Einstein condensate. 

In the presence of interactions, condensation can be 
understood as effective single-particle states (natural or- 



bitals) being macroscopically populated. Following Pen- 
rose and Onsager one needs to diagonalize the one- 
particle density matrix pij and study its eigenvalues 



(18) 



(i) If only one of them scales proportionally to the to- 
tal number of bosons in the system (iV&), one has sim- 
ple Bose-Einstein condensation, (ii) If more than one of 
them scales proportionally to Nb, this means that frag- 
mentation occurs [5ll ]. (iii) One can also have many of 
them that scale proportionally to Nff (0 < a < 1), and 
one would then say that quasicondensation occurs, (iv) 
Finally, if all of them scale as 0(N°), the system does 
not exhibit Bose-Einstein condensation or quasiconden- 
sation. Interestingly, exact results for expanding hard- 
core bosons in one dimension [l6j have already shown 
that the largest eigenvalue of the one-particle density ma- 
trix in a system out of equilibrium can diverge (ex yWfe 
in that case), while still being composed by bosons with 
many different momenta. In those systems, there were 
as many populated momenta as those available in the 
Fermi sea of the noninteracting fermions to which hard- 
core bosons were mapped. 

Within the mean-field approximation, the elements of 
the one-particle density matrix reduce to 



Pi 



(19) 



The expression in Eq. ([TO]) , which holds because of the 
intrinsic form of the Gutzwiller state in Eq. implies 
that the one-particle density matrix can only attain at 
most one eigenvalue which grows as the total number of 
particles grows; thus if fragmentation or quasicondensa- 
tion occurs in a real system, the Gutzwiller mean-field 
approach will fail to reproduce them [52| . A scenario in 
which all A's are 0(N®), on the other hand, can of course 
be reproduced by this approximation (e.g., the Mott in- 
sulator at zero temperature and the normal phase at high 
temperature). 

In Fig. O we show the time evolution of the largest 
eigenvalue Ao of pij for three different values of U/J 
{U/J = 25, 50, and 100) in our systems. Due to the 
presence of the superfluid rings in the initial state, Ao is 
not always small at t = 0. (As expected, its value at 
t = decreases as U/J is increased.) For all the sys- 
tems studied in this work, we find that the occupation of 
the highest populated natural orbital always grows dur- 
ing the expansion (Fig. \§§ , while the occupation of the 
other natural orbitals remains negligible at all times. In 
addition, we have verified that Ao always increases pro- 
portionally to the total number of particles. Hence, at 
least within the mean-field approximation, we find that 
one natural orbital is macroscopically populated. This 
means that one of the scenarios (i)-(iii) could be relevant 
to the physics of the real systems. Experiments and/or 
other theoretical approaches will need to be used to iden- 
tify which of those scenarios (if any) is the correct one. 
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FIG. 9. (Color online) Comparison between the time evo- 
lution of the occupation of the highest populated natural 
orbital Ao for systems with U/J = 25 (black), U/J = 50 
(red), and U/J = 100 (green), after turning off traps with 
V/J = 0.053, 0.162, and 0.356, respectively. The time is 
reported in units of h/J. 

VI. CONCLUSIONS 

We have presented a detailed study of the expansion of 
n = 1 bosonic Mott insulators in two-dimensional optical 
lattices. In all cases studied here, the initial states were 
prepared in such a way that the largest fraction possible 
of the bosons would be part of the Mott-insulating do- 
main. Two setups were considered: the expansion along 
one direction of an anisotropic system, which was initially 
trapped only along the x direction, and the expansion of 
an isotropic system, which was allowed to expand in all 
directions. In the first case, we have shown that the ex- 
pansion of the Mott insulator leads to a high population 
of two well-defined momentum states. Those momenta 
were found to be controlled by the ratio of the hopping 
amplitudes along the x and y directions. In the second 
case, the expansion of a Mott insulator in all directions 
produced a simple condensate composed by bosons with 
many different momenta. Given the approximated na- 
ture of our calculations, experiments, and/or unbiased 



theoretical work, will need to verify our findings and 
address the nature of the condensation seen within the 
Gutzwiller mean-field approximation. For the two set 
ups mentioned before, we also presented analytical argu- 
ments to determine which momentum modes are popu- 
lated during the expansion. 

It is important to notice that the phenomena we have 
observed here arise because of a subtle interplay between 
quantum tunneling and strong correlations. Both effects 
play a crucial role in the generation of the condensates 
as, on the one hand, tunneling is necessary for the system 
to expand, and on the other hand, interactions are neces- 
sary for the bosons to redistribute the energy and couple 
coherently. In the absence of interactions, the system 
still expands, but the momentum distribution remains 
unchanged from its initial value, i.e., it remains feature- 
less if one starts from purely localized particles in the 
lattice sites (a Fock state), and no condensation occurs. 

Finally, our study supports the conclusion in Ref. [HI , 
and generalizes it to the soft-core boson case, that 
strongly correlated atom lasers can be generated from the 
expansion of Mott-insulating states. Furthermore, after 
the results in Sec. IIV1 one can see that the momenta of 
those matter waves can be controlled by changing the 
ratio between the hopping parameters in the different di- 
rections in the lattice. In addition, if one cannot change 
those, in an isotropic lattice (Sec. W\ < the momenta of 
the lasers can be controlled by changing the direction 
along which the Mott insulator is allowed to expand. In 
the latter setup, the best results would be obtained by 
constraining the Mott insulator to expand along the di- 
agonals in the lattice. 
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